%%%%%%%%%%%%%%%%%% Poisson distribution CI2 %%%%%%%%%%%%%%%%%%%%%%%% 
clear

n=1;
m=100;  
r=50;
k_item=[0:m];
lambda_range=60; 
alpha=0.05;
alpha=0.17;
beta =0.9;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CI2
p1_tilde=0.5*chi2inv(  alpha/2, 2*k_item  )/n; 
p1_tilde(1)=0;
p2_tilde=0.5*chi2inv(1-alpha/2, 2*k_item+2)/n;
p_tilde=[p1_tilde; p2_tilde];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%找出I,J=>解出P%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

II=[];
JJ=[];
P=[];
for k=0:r
    
     p1=p_tilde(1,k+1);
     p2=p_tilde(2,k+1);
     
     H1=poisscdf(m,p1)-poisscdf((k_item-1),p1)-(1+beta)/2;
     %H1(1)=poisscdf(m,p1)-(1+beta)/2;
     H2=poisscdf(k_item,p2)-(1+beta)/2; 
     
     c1=min(H1(H1>=0));
     c2=min(H2(H2>=0));

     %if k==0; 
     %    I=0; J=0; 
     %else 
         I=find(H1==c1)-1;
         I=max(I);
         J=find(H2==c2)-1;
         J=min(J);

  %   end
     II=[II,I];
     JJ=[JJ,J];
   
     poiss=[];
     %for x=(I+1):J
     for x=I:J
         strf=[ 'exp(-v)*v^' num2str(x) '/'  num2str(gamma(x+1)) '+'];
         poiss=[poiss,strf];
     end
     poiss(end)=[];
     poiss=[poiss '-' num2str(beta)];
%----------------------     
     fun = fcnchk(poiss);
     a=feval(fun,0);
     b=feval(fun,lambda_range);
     poiss(poiss=='v')='x';
     
     
     if(a*b<0) %一個解
          p=fzero(poiss,[0 lambda_range]);
          p=[p NaN];
     else     %無解或兩個解
          for i=0:1:lambda_range,
              if(feval(fun,i)>0)
                  break;
              end
          end
          if(i==lambda_range)|(i==0)
              p=[NaN NaN];
          else
             p1 = fzero(poiss,[0 i]);
             p2 = fzero(poiss,[i lambda_range]);
             p=[p1 p2];
          end
     end 
     P=[P;p];
end
sprintf(' m     k        I        J       v') 
%output= [ m*ones(r,1) (1:r)' II' JJ' P]
output= [ m*ones(r+1,1) (0:r)' II' JJ' P]
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%對解出的P做排序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:length(output(:,1))
   if isnan(output(i,6))==0
      output=[output; output(i,1:4), output(i,6), NaN]; 
   end
end
[B,IX]=sort(output(:,5));
sprintf('        n         k        I        J       Solve P   ' ) 
output=output(IX,1:5)
%output=[n 0 0 0 NaN; output];
%output=[n 0 0 0 NaN; output; n n n n NaN];

%%%%%%%%%%%%%%%%%對某個P..找出對應的x, 算出Conf._Level.%%%%%%%%%%%%%%%%%%%%%%%% 
aa=[];
for i=1:length(output(:,1)) %  i=2 ;
   lambda=output(i,5);
   x=[];
   for j=1:length(output(:,1)) % j=4;
        k=output(j,2);  I=output(j,3); J=output(j,4);
        
      if(I==0)
        temp=sum(poisscdf(J,lambda));
      else 
        temp=sum(poisscdf(J,lambda))-sum(poisscdf(I-1,lambda));
        %temp=sum(poisscdf(J,lambda))-sum(poisscdf(I,lambda));
        %temp=poisscdf(J,lambda)-poisscdf(I,lambda);
      end           
        
      if  temp-beta>0   
          x=[x;k];
      end
   end    
   a=sum(poisspdf(min(x):max(x),lambda));
   aa=[aa,a];  
end
minconf=min(aa)
%minconf=min(aa(2:length(aa))); % 排除k=0 沒有解的情況
sprintf('      n  observation  tolerance  tolerance   theta    coverage \n                        lower      upper             probability   \n                        bound      bound  ') 
output=[output,aa']


%%%%%%%%%%%%%%%%%%%%%%%%%計算ave. cov.prob.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

P=[];
bd=[];
tempbd=[0;output(:,5);lambda_range];
%tempbd=[0;output(:,5);1];------------------------
for i=1:length(tempbd) 
   if isnan(tempbd(i))==0
      bd=[bd; tempbd(i)]; 
   end
end

for i=(1:length(bd)-1)
   tempP=(bd(i)+bd(i+1))/2;
   P=[P;tempP];
end

aa=[];
result=[];
for i=1:length(P) %  i=12;
   lambda=P(i);

   x=[];
   for j=1:length(output(:,1)) % j=1;
        k=output(j,2);  I=output(j,3); J=output(j,4);
        
      if(I==0)
        temp=sum(poisscdf(J,lambda));
      else 
        temp=sum(poisscdf(J,lambda))-sum(poisscdf(I-1,lambda));
        %temp=sum(poisscdf(J,lambda))-sum(poisscdf(I,lambda));
        %temp=poisscdf(J,lambda)-poisscdf(I,lambda);
      end  

        if  temp-beta>0;    
            x=[x;k];
        end
   end    
   a=[lambda,min(x),max(x)];
  
   
   if length(a)==1
       aa=[aa;a NaN NaN];
       resultt=0;
   else 
       aa=[aa;a];
       poiss=[];
       for x=a(2):a(3)
              strf=[ 'exp(-v)*v^' num2str(x) '/'  num2str(gamma(x+1)) '+'];
              poiss=[poiss,strf];
       end
       poiss(end)=[]; 
% average coverage probability for bounded parameter space (bound1, bound2)
       bound2=9;
       bound1=0;
       newdl=max(bound1,min(bd(i),bound2));
       newdu=max(bound1, min(bd(i+1),bound2));
       resultt=1/(bound2-bound1)*int(poiss,newdl,newdu);
%resultt=1/(bound2-bound1)*int(poiss,newdl,newdu); 
      % resultt=int( poiss,bd(i),bd(i+1));
   end   
   
   result=[result,resultt];    
   
end
sprintf('   inter.p     min(x)      max(x)   ' )
aa 

%average coverage probability  for parameter space (0, lambda_range) 
%average_coverage_probability=(double(sum(result))/lambda_range)

%average coverage probability for parameter space (0, lambda_range) 
average_coverage_probability=(double(sum(result)))

result=sprintf('m= %g,   r=%g,   minimum coverage probability =%.4f,  average coverage probability= %.4f',m,r,minconf,average_coverage_probability)
% result=[m r minconf double(sum(result))/lambda_range]